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Fluctuating hydrodynamics for dilute granular gases: a Monte Carlo study 
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We investigate hydro-dynamic noise in a dilute granular gas during the homogeneous cooling 
state, by means of a proper application of the Direct Simulation Monte Carlo (DSMC) algorithm. 
The DSMC includes a source of randomization which is not present in Molecular Dynamics (MD) 
for inelastic hard disks. Notwithstanding this difference, a fair quantitative agreement is found, 
including a violation of the fluctuation-dissipation relation for the noise amplitude of the same 
' f—Z • order observed in MD. This study suggests that deterministic collision dynamics is not an essential 

ingredient to reproduce, up to a good degree of approximation, hydrodynamic fluctuations in dilute 
£— n , granular gases. 

<& : 

^ ; I. INTRODUCTION 

(n: 

Granular gases represent an important benchmark for theories and methods of non-equilibrium statistical mechan- 
ics [l(. A granular gas is the idealized counterpart of fluidized granular media, i.e. collections of N ~ 10 2 -=- 10 3 
macroscopic grains (typically spheres with diameters a in the range 10 -j- 10 -3 m), contained in volumes V such 
that the packing fraction, or occupied volume, is smaller (often much smaller) than 10%: note that the small packing 
fraction ~ Na d /V is still compatible with a non-negligible collision rate, ~ Na^ 1 /V , where d is the dimensionality 
of the system. Within these conditions, the main difference with respect to an ideal gas is given by inelastic collisions: 
kinetic energy is lost during an (hard core) interaction, ruling out many basic assumptions of equilibrium statistical 
mechanics and leading to several phenomena such as: breakdown of energy equipartition [2-4], spontaneous symmetry 
breaking (vortices and clustering) [5|, Q , "Maxwell demon" [7| , "ratcheting" effects [8[ , and lot more. In both theo- 
retical investigation and in real world applications, granular "particles" are often found in a small number, thousands 
is a common order of magnitude, as mentioned before. This makes fluctuations a fundamental ingredient of granular 
f_\ ' statistical descriptions 6]. While continuum theories, e.g. granular hydrodynamics and kinetic theory, describing local 
averages for density, velocity and energy and their evolution, have achieved an advanced degree of maturity Q, a 
good description of noise accompanying such averages is still an open problem. Near equilibrium, Einstein fluctuation 
theory, resumed by Landau and Lifshitz in their fluctuating hydrodynamics [lOj . based on equilibrium fluctuation- 
dissipation relations (FDR) [IJ, [l2|, simplifies the description of noise for transport processes. In granular gases, 
fluctuations do not satisfy, in general, FDR. The simultaneous presence of external random drivings and diluteness, 

t— I • usually guarantees the validity of FDR [13, [Ij]. FDR is, instead, violated at packing fractions > 10% [15|, Qj| and in 
a non-driven granular gas, even dilute, in the so-called Homogeneous Cooling State (HCS) [17L Il8j| . 

iy-N . A recent study on dilute inelastic hard particles in the HCS, with Molecular Dynamics (MD) simulations [l8[ and 

then an analytical treatment [19| , have shown that fluctuations of the hydrodynamic transverse velocity field deviate 
from the classical Landau-Lifshitz theory for two reasons: 1) noise is not white, 2) it does not satisfy the FDR of 
the 2nd kind. The analytical study is based on projection operator formalism, starting from the "pseudo-Liouville" 
operator for the deterministic evolution which only includes free flights and instantaneous hard-core inelastic collisions 
at contact. The only statistical averaging is made with respect to an ensemble of initial conditions, while several 
hypothesis are introduced: among the others, spatial homogeneity and a special time-scaling (inherent to HCS), for 



one-particle and two-particle distributions, such that all time-dependence is contained in the "granular temperature" 
T g = m(v 2 )/d. The last part of this analysis makes use of the Molecular Chaos assumption, which is reasonable for 
dilute systems. Anyway, one wonders which ingredients of this sophisticated analysis are essential and which can be 
neglected. For instance, the use of Molecular Dynamics (MD) for simulations and of the rigorous manipulation of the 
pseudo-Liouville operator as starting point of the analytical treatment let us think that the deterministic nature of 
the hard spheres dynamics is crucial. 

Here we present a different model, which in principle could produce different results, but in the end reveals to be 
in fair agreement with the previously cited study. It is a stochastic model, where collisions are treated randomly, 
in sharp contrast with deterministic dynamics employed in Molecular Dynamics simulations and described by the 
pseudo-Liouville operator. It consists of the classical "Direct Simulation Monte Carlo" (DSMC) algorithm [2fJ in 
its non-homogeneous form, i.e. placed on a spatial grid to measure modes of the transverse velocity field, with the 
use of "Active particles" to guarantee Molecular Chaos (the algorithm is explained in the Appendix). The inelastic 
collision rule, relating post-collisional velocities to pre-collisional ones, and the time-rescaling are the only ingredients 
which are conserved with respect to the MD model. Our question is: are those ingredients sufficient to reproduce the 



hydrodynamic noise properties of the deterministic model? 

Similar studies have been performed for hydrodynamic fluctuations in the elastic case [2l|, [22[ . More recently, the 
DSMC approach has also been applied to fluctuations of the global energy fluctuations (a homogeneous, not spatially 
dependent quantity) in inelastic systems [23, 24], obtaining a good agreement with the amplitude of fluctuations 
measured in MD simulations and explained in terms of projection operators applied to the deterministic dynamics |25| . 
A similar conclusion was available for that study: stochastic "Monte Carlo" treatment of collisions is sufficient to 
reproduce MD noise properties. A successive study |26| also showed how certain quantities (cumulants of fluctuations 
in the stationary state) may be obtained without resorting to projection techniques. 

In the present study, in a sense, the observation is even more striking: the exchange of particles among differ- 
ent "copies" , needed to reduce the incidence of re-collisions and finite size correlation effects, is a strong source of 
randomization. Surprisingly it does not affect the hydrodynamic fluctuations. 

The organization of the paper is the following: in Section[H]wc first describe the model, then in Section|TTT]we present 
the results, first verifying Molecular Chaos and then measuring the hydrodynamic noise in time-autocorrelation and 
amplitude. Finally we draw Conclusions. We have dedicated an Appendix to a detailed description of the DSMC 
algorithm used for this work. 

II. THE THEORY 

A. Definition of the model 

This Direct Simulation Monte-Carlo study approximates the dynamics of a dilute system of N smooth inelastic 
hard disks of mass m = 1 and diameter a. The inelastic collisions between the disks i and j change the particles 
velocities following the relations 






1 + a 
v i = Vj + — — {<r-Vij)a (1) 



where v' corresponds to post-collisional velocity, a is the coefficient of restitution and Vjj = Vj — Vj is the relative 
velocity. The term j is a random unit vector (see Appendix). The size of the system is L x L and the boundary 
conditions are periodic. 

The main hypothesis used here is that the system, starting in a uniform equilibrium state with granular temperature 
T g (Q) — To, evolves to the Homogeneous Cooling State (HCS), which is characterized by a single time-scale measured 
by the temperature T g (t): any other quantity depends on time only through T g (t). Apparently, as observed in many 
previous studies [23], in a homogeneous setting, Molecular Chaos is sufficient to guarantee this hypothesis. 

The HCS is unstable against spatial fluctuations: this instability appears at scales larger than a critical length L c (a), 
therefore it can be avoided by taking the linear size of the system L < L c [28] . It is possible to analyze the effects 
of spatial fluctuations by deriving mesoscopic equations through a linearization around the HCS |5j. The resulting 
equations are generally coupled, but in the Fourier representation the transverse velocity field results decoupled 
from the other modes. We are interested, in particular, in the fluctuations of its largest mode, i.e. of the smallest 
wave-number |k m | = 2ir/L. We choose the wave vector k m parallel to the axis x and the quantity above results to be 

U ± (t) = f2v y>j (t)e^(i^^) (2) 

j 

where v y j(t) is the component of the velocity of the particle j in the direction y and Xj(t) is its coordinate along the 
x axis. We want to verify the hypothesis that the fluctuations of U±(t) obey a linear Langevin equation, [5l Il8j 



dtUxit) = -v{t)k 2 m U±{i) + y/vth(t)Nc v(t)kU(t) (3) 

where v(t) is the kinematic viscosity which, in the HCS, is proportional to \JTg (see |9( for definitions), and Vth{t) = 

The last term in the Eq. Q describes an internal noise corresponding to the rapid (microscopic) degrees of freedom 
of the system. This complex noise is assumed here to be white, Gaussian and with correlations given by 

me(t r )) = 6{t-f) (4) 

WW)) - (5) 



with £*(£) the complex conjugate of £(£). Deviations from the white noise assumption are expected [lj|, but they are 
in general quite small. The most evident consequence of inelasticity is, instead, the variation of the dimensionless 
coefficient c in Eq. ((3]), which is equal to 1 at equilibrium (a = 1); its departure from c = 1 represents the fact that 
the fluctuations of U± (t) do not satisfy the FDR. 

Since the system is cooling, the typical velocity of the particles becomes smaller and smaller, leading to increasing 
rounding errors and other numerical problems. These problems can be avoided with the use of a procedure which, in 
the HCS, maps the dynamics to a steady state by means of a time-rescaling [23, (29|. Such procedure is explained in 
the next subsection. 

B. The stationary representation 

In the HCS, the granular temperature obeys the following equation: 

d t T g (t) + t H (t)T g (t)=Q (6) 



where £# is the cooling rate that results proportional to y/T g (t). 

The stationary representation of the HCS |27| consists in introducing a new time scale r defined by 

loot = In — (7) 

to 

with ujq ant to arbitrary constants, implying the definition of rescaled velocities v(r) = v(t)u)ot. It is easy to see [27] 
that observing the system on this new time scale is equivalent to apply a positive continuous drag to all particles 
9 t v(t) = wqv(t). This naturally leads to define also the rescaled analogous of U±(t): W±(t) = U±(t)uot. 
We tune wo in order to have T g (r —^ oo) — Tq (so that the length of transients is reduced), a result which is obtained 
by taking ujo = £#(0)/2. In the steady state, we simplify the notation using vth — uth(O). With these choices, the 
Langevin equation (j3|) is mapped onto a new equation 



d T W ± (r) = -^kl-^W^r) + ]JNvf h ^ £(r), (8) 

where Ao = L 2 /(No) = l/(na) is proportional to the mean free path and v = u(t)/[XoVth(t)] and k m = k m Xo are 
dimensionless rescaled viscosity and wave number, respectively. 

If vk^ n — luq < 0, which is equivalent to the condition of stability of shear modes in the HCS, the above equation 
leads to a dimensionless autocorrelation function given by 

Cl|T)5 wp4( lt 4-./. (9) 

v th 2 V J 

with Tq = vk^yth I 'Ao — ^>o the characteristic time of decay. 

Based on Eq. ©, we can obtain the dimensionless kinematic viscosity v and the "FDR violation" coefficient c from 

a measure of C±(t): 

v th k 2 ^T ) 

2CM0) 



iv(T u; + l) 
i.e. v and c are obtained measuring the amplitude Cj_(0) and the decay time tq 



III. NUMERICAL RESULTS 



The aim of this section is to measure v and c by means of Equations (flQ)) and (TiT]) in DSMC simulations, comparing 
them to the results of MD. For this purpose a DSMC algorithm must be carefully devised to guarantee the validity 
of the HCS assumptions and to allow for the measure of W±(r). This is done in the next subsection. 



A. Stationary representation for the inhomogeneous DSMC 

The DSMC algorithm used here is explained in details in the Appendix. The essence of a DSMC algorithm, which 
has been originally designed as a tool to numerically solve Boltzmann equations, is a stochastic computation of 
collisions: pairs of particles and their orientation vector a are chosen randomly, with probabilities dictated by the 
Boltzmann collision integral. In homogeneous configurations, a pair of particles is chosen among all N particles of the 
system, disregarding spatial coordinates: the large number N guarantees Molecular Chaos. Here, to avoid shear mode 
instabilities, we keep the total size of the system under the critical size L c , expecting to observe a homogeneous regime. 
Nevertheless, the quantity under scrutiny, U±(t), requires that collisions are treated with a good spatial resolution, 
i.e. colliding particles must be close to each other to give the correct contribution to the variation of U±(t). For 
this reason, a inhomogeneous DSMC algorithm is necessary. The system is partitioned in m c non-overlapping cells 
of size l c : during the free-streaming step each particle can move from a cell to a neighboring one, but during the 
collision step the particle can collide with particles in the same cell only. The resolution for the measure of U± (t) is 
improved increasing m c ; at the same time, when m c increases, the average number of particles in each cell N c = N/m c 
decreases, threatening the Molecular Chaos assumption. In principle the perfect resolution would be achieved with 
l c ~ er, but this would result in a number of particles per cell N c ~ na 2 which, for diluteness, is required to be much 
smaller than 1. We will see in the following (Fig. [2]) that taking l c ~ 15er is sufficient to restore a good resolution for 
measuring U±(t). Such a choice, however, at the chosen packing fraction na 2 = 10~ 2 gives N c ~ 2.5, which is so low 
that the Molecular Chaos could be invalidated. A first check is to study the sensitivity of the rescaling procedure © 
(tested in [27[ in the homogeneous DSMC simulations and in [30]) to the parameter m c (or equivalently N c ). 

Having this aim, we have characterized the entire system analyzing its granular temperature T g as function of the 
rescaled time r for different choices of the cell number m c at a and N fixed (see Fig. [1]), comparing it with MD 
results. The time-rescaling is applied with a choice of wo such that, in the steady state it is expected T g =T . From 
the simulation data it results that a steady state is always reached, i.e. the granular temperature is constant in all 
cases. A homogeneous DSMC (i.e. m c = 1) gives an agreement with the expected value which is better, of a few 
percentages, with respect to MD. When m c is increased, however, the DSMC gives very bad results. This is due to 
the small value of N c , which results in strong finite size effects, for instance fake re-collisions which invalidate the 
Molecular Chaos hypothesis. 
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FIG. 1: The rescaled granular temperature T g as function of the rescaled time, in DSMC simulations with different / factors 
and cell numbers m c , compared with MD simulations, in the case of JV = 1000, a = 0.8 and na 2 — 0.01. 



The way of solving this problem is to simulate a larger number of particles N = fN with / > 1, keeping the 
collision statistics per particle as that of a A-particles system (see the Appendix for a precise description). This is 
equivalent to simulate / copies of the original ./V-particles system and let particles of different copies swap at each 
time step. The effect of / > 1 is immediately seen in Fig. [TJ the steady state granular temperature perfectly agrees 
with that measured in MD, as soon as / ~ 10, even if m c = 225. Much larger values of / appear to slightly improve 
the value of T g . On the other side a too much large value of / can be very expensive in terms of cpu-timc. 

A crucial quantity, here, is the number of virtual particles per cell N c = fN c = fN/m c . If m c is increased (for 
instance in Fig. Q] it is changed from 225 to 400), / should be increased. Apparently it is sufficient to have N c ~ 25 
to have a good agreement, for T„, between DSMC and MD. 
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FIG. 2: The rescaled correlation function C±(t) (see Eq. [9} as function of the rescaled time, obtained from DSMC simulations 
with different / factors and cell numbers m c , and from MD simulations, in the case of N — 1000, a = 0.8 and na = 0.01. 



As discussed above, the necessity of a large number of cells comes from the sensitivity of U±(t) to the spatial 
localization of collisions. This is shown in Fig. [2J where the decay of C±(t) is displayed: the homogeneous DSMC 
(m c = 1) shows an exponential decay with a t completely different from the one observed in the MD simulations. 
This behavior is insensitive to changes of / from 1 to 10. In order to obtain the correct decay, it is necessary to use 
the inhomogeneous DSMC approach with m c 3> 1. A partition of m c — 20 x 20 for N — 1000, in fact, reveals an 
exponential decay very close to the MD results: also in this case, a change of / from 10 to 100 is unimportant. 
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FIG. 3: The rescaled correlation function C±(t) as function of the rescaled time, from DSMC simulations with different total 
particle numbers iV and same average number of particles per cell -/V c . Data from MD simulations are shown for comparison. 
The parameters are the same as in Fig. [TJ 



When N is increased at fixed density (and mean free path), the size of cells must be kept constant, which implies 
m c oc N . If / is also kept constant, N c is automatically preserved and we expect to have a good comparison with MD 
simulations in terms of Molecular Chaos, as well as decay of C±(t). This situation is fairly verified in Fig. O where 
the decay of Cj_(t) is shown to be always close to the MD results. 

This approach is valid in general and, with proper values of / and m c , we can compare, for different values of a 
and N, the properties of the fluctuations of W± between DSMC and MD simulations. Guided by this preliminary 
analysis, we have set the value of / = 10 and N c — 25 for all the following discussion. 



B. Analysis of the fluctuations 

In the previous section we made the hypothesis that the fluctuations of the slowest mode of the transverse velocity 
field U±(i) follow the Langevin equation ([3]), or ((5J) in the rescaled representation: this implies that its correlation 
function follows the exponential decay in Eq. ©, with characteristic time Tq. 
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FIG. 4: The dimensionless kinematic viscosity v, obtained from Eq. (|10|l and from the exponential decay of the correlation 
function in Eq. JSJ, as function of the inverse of the number of particles N for different values of a. The circles and squares 
correspond to data obtained from DSMC an MD simulations, respectively. The solid and dashed lines are the best linear fit. 
The density is the same as in Fig. [T] 



In order to apply the correct thermostat v = ojq'v, and to obtain the kinematic viscosity v from Eq.©, it is 
necessary to know the value of u>q = £#(0)/2. This is computed, in the homogeneous case, giving the theoretical 
expressions [27| 



Cj?(o) = 



vth(o) n? 



l + —a 2 {a) 



where the coefficient 02(0;) is 



02(a) 



16(l-a)(l-2a 2 ) 
57 - 25a + 30a 2 - 30a 3 



(12) 



(13) 



We have verified that the prediction in Eq. (|I2|) is good also for the inhomogeneous DSMC and for MD. 

Eq. © is expected to be valid in the large N limit: for this reason we have performed simulations for different 
values of N (1000, 1500, 2000, 2500 and 3000), fixing the particle density n of the system, and measuring the time 
decay to. The data obtained using (TTOl) and (fl~2"|) in DSMC and MD, are shown in Fig. [4] The dimensionless viscosity 
v appears to be a linear function of N , therefore a best linear fit allows to extrapolate the asymptotic value for 
N — >• 00. The agreement between MD and DSMC data is good in all cases, even at not too high values of N. 

Our second question is if the inhomogeneous DSMC can reproduce the MD results for the noise term in the Langevin 
equation ([3]) or ((SJ, in particular the quantity v' = cv that is the viscosity term appearing in the noise. Analogously 
to Fig. IU we have analyzed its dependence versus iV^ 1 . Also in this case we have found a linear law, shown in 
Fig. [5] a small discrepancy is observed, for this quantity, between DSMC and MD data, when N is increased. The 
disagreement results more evident if the collisions are more inelastic. Our interpretation of this small difference is 
that the DSMC dynamics at small time-scales is slightly different from MD: indeed, this disagreement is less evident 
at large time-scales, for instance in the measure of v. If we consider also that the DSMC and MD results about the 
correlation of the energy fluctuations are very similar [23| , we tend to conclude that this small discrepancy could be 
due to the use of small cells (large m c or equivalently small N c ), which is necessary in the measures of Cj_(0). A more 
general discussion of the cause of this discrepancy is given at the end of this section. 

Notwithstanding this slight discrepancy, the value of c = v' fv results significantly larger than 1 also in DSMC. The 
inhomogeneous DSMC is then able to reproduce the breakdown of the FDR as already found with the MD approach 
[18|. This is evident from Fig. [51 where the v and v' are plotted as function of a for different values or the particle 
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FIG. 5: The dimensionless viscosity v — cv obtained from Eq. (|11[) and (|10|l (see text), as function of the inverse of the 
number of particles N for the following values of a: 0.8, 0.85, 0.9 and 0.95. The circles and squares correspond to data obtained 
from DSMC an MD simulations, respectively. The solid and dashed lines are the best linear fit. The density is the same as in 
Fig.!] 



number and for TV — > oo. The extrapolated results for N — > oo, for v are in good agreement with the MD data and 
with the theoretical analysis of [19| . while those for v' still underestimate them, as shown in the fourth panel of Fig. 

m 

The ratio c is shown in Fig. [7J it clearly displays strong deviations from 1 (violations of FDR), which reach 12% 
at a = 0.8 also in DSMC. The dot-dashed curve represents the result for c from the theoretical analysis of [19| with 
the additional assumption of white noise. We note a close comparison with the data from DSMC, suggesting that 
the effect of this protocol is an effective time-decorrelation of noise. Such a conclusion should be taken with care, in 
view of other simulations, performed with different choices of / and m c : these results are not shown here, since have 
a quite high statistical uncertainty due to the large computational time required (increasing / and m c results in an 
increase of fictive particles N). They suggest that, when / and m c are increased keeping N c = fN c constant, the 
value of c stays substantially unchanged. On the other side if N c is increased (e.g. by increasing / at constant m c ) the 
value of c slightly increases and seems to slowly tend toward the corresponding MD value, which is also compatible 
with the theoretical analysis with a colored noise. The presence of time-correlation in the noise, even in DSMC, is in 
agreement with the results for C±(t) shown in Fig. [2] and [3] which are not straight exponentials, but some bending 
can be observed at small r. 

Apart from these small discrepancies, the main fact is that the inhomogeneous DSMC is able to reproduce the 
violation of FDR, which was not obvious. In fact, it is reasonable and well verified that average values (e.g. transport 
coefficients) obtained in DSMC agree with dilute MD simulations; anyway it is less trivial to observe good agreement 
for fluctuations, which - in non-equilibrium situations - can be much more sensitive to the detailed mechanisms of the 
dynamics. 



IV. CONCLUSIONS 



In this work we have obtained two main results. First, we have individuated which are the basic ingredients - in 
a Direct Simulation Monte Carlo algorithm - to correctly measure the time-decay of the autocorrelation function of 
fluctuations of a hydrodynamic field, keeping valid the assumption of Molecular Chaos: the volume must be divided 
in small cells for the purpose of collision computation, even if the system is assumed to be spatially homogeneous; 
the good resolution to appreciate the "real" time-decay of C± (r) is achieved only when N c ~ 2 which is too small to 
avoid fake recollisions: for this reason, a number (order ~ 10) of "virtual" copies of the system is necessary to restore 
Molecular Chaos. This procedure is similar to that adopted in [31(. The study of the transverse velocity field could 
also be done by dividing the system in slides, as in [2ll.l24j|. Such a procedure can be used to analyze U±(t), while we 
have no knowledge of its adaptation to other space-dependent fields. Note that - when DSMC is used to study the 
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FIG. 6: The dimensionless viscosity v — cv (empty symbols) and the dimensionless shear viscosity v (full symbols) as function 
of the coefficient of restitution a for different values of N: 1000, 2000 and 3000. The circles and squares correspond to DSMC 
and MD data, respectively. The data for N — ► oo are extrapolated, with the respective errors, from the linear fits of Fig. [3] 
and Fig. [5] The solid and dashed lines are the theoretical predictions [19| . 
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FIG. 7: The ratio c = v' /v as function of the coefficient of restitution a for N — > co (see fourth panel of Fig. [6]). The circles 
and squares correspond to DSMC and MD data, respectively. The solid is the theoretical prediction (191 ]. The value of reference 
c = 1 (dashed line) is obtained if FDR are satisfied. The dot-dashed curve is the result of the theoretical analysis of [ly] with 
the additional assumption of white noise. 



homogeneous cooling state, i.e. keeping the system size below the critical size for instability, L < L c - one usually 
does not divide the system in cells, and still obtains excellent agreement for one-point observables. Previous study of 
cooling granular gases with inhomogeneous DSMC [30L l31j where in fact proposed to analyze the departure from the 
HCS. 

Our second result concerns a quantitative agreement with Molecular Dynamics results, which is very good for 
the exponential tail of the autocorrelation of fluctuations. The amplitude of fluctuations, on the other side, shows 
small discrepancies with respect to MD results, which become more evident as a is reduced and iV increases. In 
the extrapolated N — > oo limit, these discrepancies can reach (for a — 0.8) an underestimation of order 10%: 
nevertheless it is still possible to appreciate a ~ 12% of deviation from the validity of FDR, which is one of the most 
peculiar property of granular hydrodynamic fluctuations in the HCS. In DSMC simulations, deviations from a perfect 
exponential decay of C±(t) are also visible, at small times, but these are only qualitatively similar to those observed 



and predicted in [19(. 

Hydrodynamic noise is always important in granular systems, which - in terms of number of elementary constituents 

- are much smaller than molecular fluids (10 4 instead of 10 20 particles). Future research is needed to understand the 
properties of granular hydrodynamic noise in steady state models: even if we have shown that these properties do not 
depend on the detailed collisional mechanism, it seems that they depend on the choice of energy driving protocol [32| . 
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Appendix 

In this Appendix we give a detailed description of the DSMC algorithm used in this work. The algorithm involves 
the simulation of N = fN Active particles with / > 1. This is a trick often used in DSMC: for molecular gases - where 
N ~ 10 20 - it is customary to use / <C 1 in order to have a manageable number of Active particles. In granular gases, 
where the number of real particles is small, one uses / > 1 in order to have enough particles in a cell to guarantee 
Molecular Chaos. This is our case. 

As initial condition we have chosen a uniformly random spatial configuration of the particles and a Gaussian 
distribution of their velocities. After initialization, the dynamics consists of a cycle with the following two steps: 

• Streaming step (evolution of positions and velocities ignoring possible collisions) 

• Collision step 

The streaming step consists in computing, for every particle i, the time-discretized version (with a fixed At, small 
enough) of the evolution equations 

r* = v, (14) 

Vi = uj Vi (15) 

For the purpose of the collision step, the system is partitioned in m c equal non-overlapping cells of size l c (an 
optimal value l c ~ 15<r has been individuated, as discussed in Section IlII A[) . The average number of fictive particles 
in the cell i is then N c = N/m c . If m c is too large, the number of real particles in a cell N/m c , results too small and 
invalidates the Molecular Chaos hypothesis, due to fake recollisions. This is avoided by taking / > 1 large enough (a 
value of ~ 10 seems optimal). 

In each cell k the local average collision frequency - per particle - Uk = &rik y/^Tk/m is calculated, based on the real 
density n,k = Nk/{fl1) where Nk is the number of fictive particles in the cell and an estimate of the local temperature 
Tfc (i.e. the variance of the local particle velocity distribution). Then a number of collisions NkUJkAt is performed: this 
is obtained choosing, at random, pairs of particles i and j, and unit vector (uniformly distributed angle) a, accepting 
them with a probability proportional to (vi — Vj) • a. The accepted pair is updated with equations ([1]). 

The algorithm described here is the classical DSMC algorithm as proposed by Bird [20j (variant exist which include 
fluctuations of the number of collisions with the introduction of an internal clock for each cell) . It is interesting, in 
this particular case, to realize that the number /, when integer and larger than 1, can be thought as a number of 
"copies" of the system. In fact, if At is small enough to guarantee that NkOJkAt <C Nk for any cell k (which is always 
true in our simulations), then one can always virtually separate the Nk fictive particles of the cell into / groups of 

- averagely - Nk particles, such that collisions occur only between particles of the same group. Given that colliding 
pairs are chosen at random among the Nk particles, at each step this virtual separation into / groups is done in a 
new way, which is equivalent to say that - at each step - some of the particles move from one group to the other, 
remaining of course in the same cell. 
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This is just an alternative representation of the algorithm, which underlines the randomizing mechanism - intrinsic 
of DSMC and not present in MD. 
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